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Resume 



Ce memoire resume quelques concepts importants en cosmologie et presente l'etude 
faite par l'auteur sur l'origine de la lumiere intra-amas. 

Pour la realisation de ce projet, l'auteur a tout d'abord recherche dans la litterature 
les parametres a utiliser pour des simulations en langage FORTRAN dont les algo- 
rithmes de base sont, dans la premiere partie du projet, particule-particule et, dans la 
seconde, particule-particule/particule-maille. L'auteur a egalement modifie des codes 
IDL et UNIX. Enfin, le projet necessita des centaines de simulations d'amas isoles dont 
les resultats ont ete analyses en collaboration avec les membres du groupe de recherche 
et soumis pour publication (Barai, Brito & Martel 2009). 

Les resultats principaux des simulations decrites dans ce document sont: 1) la des- 
truction des galaxies naines par des fusions domine sur la destruction par des marees, 
et 2) la destruction des galaxies par des marees est suffisante pour expliquer la lumiere 
intra-amas observee. 

Finalement, les resultats d'amas isoles ont ete generalises a une region significa- 
tive de l'Univers. Ainsi, l'auteur a contribue a la mise en oeuvre d'une simulation 
particule-particule/particule-maille et a l'analyse commune des resultats obtenus a ce 
jour. Les resultats reproduisent la fonction de luminosite de Schechter, et suggerent 
que l'approche utilisee est valide et que les resultats sont robustes. 



Abstract 



This thesis presents a review of related important concepts in cosmology followed by 
details of the author's role in a research project on the origin of intracluster light. 

The author's role in the development of the simulations varied from searching pa- 
rameters in the literature, through writing and modifying code in IDL, FORTRAN, 
and UNIX to carrying out hundreds of simulations using the particle-particle algorithm 
described in this thesis, as well as partaking in joint analysis of the simulation results. 
Part of this work in the isolated cluster simulations has been submitted for publication 
(Barai, Brito & Martel 2009). 

The main results of the simulations described in this thesis are: 1) destruction of 
dwarf galaxies by mergers dominates destruction by tides, and 2) destruction of galaxies 
by tides is sufficient to explain the observed intracluster light. These results support 
the accepted explanation for the origin of the intracluster stars. 

In an ongoing, second stage of the simulation, which extends the isolated cluster 
results to a cosmologically significant region of the Universe, the author similarly as- 
sists in the implementation of a particle-particle/particle-mesh simulation and the joint 
analysis of the results to-date. The results are as per the Schechter luminosity function, 
and suggest the approach used is valid and the results obtained robust. 



Avant-propos 



The list of people to thank is endless, as usual. Thanks to everyone in my family, to be 
sure, for their patience, but within the university, I thank particularly Laurent Drissen 
for encouraging me to proceed towards the master's degree when he did — surprising 
what a word at the right time can do! Thanks to Serge Pineault for his understanding 
and encouragement along the way, especially when I had recently returned to the study 
of science. Thanks to other graduate students, such as Simon Cantin, Jean-Frangois 
Robitaille, and Veronique Petit, for answering my many questions, sometimes repeat- 
edly! Thanks to my friend Andrea Clark for copy editing so many mistakes out of this 
text — at such short notice! Thanks last (certainly not least!) to Paramita Barai and 
my Director Hugo Martel for their guidance, patience, and understanding throughout 
my time studying with them. While the above individuals (and many others) have 
helped a great deal, any error surviving to the printing stage is exclusively my own. 



To my mom, Wilma Webster, and my dad, Urbano 
Brito, who always encourage me to continue learning; 
to uncle Jim Webster, and to my siblings Cendy, Moy, 
and Steven who are also ever in my heart and mind. 



vi 



Release me for a little while 

And let me go 

Out in the twilight 

Where over red fires of sunset 

The thin moon gleams 

There let me pass 

Let mortal shadows drop away 

And leave the spinning earth 

Behind me 

In the weightlessness of outer space 
I would adventure 
Towards far galaxies 

— Andromeda Webster 
Flashes from space. 



Lo ! That some we loved, the loveliest 
and the best of all time's vintage pressed, 
have drunk their cup a round or two before 
and crept silently to rest! 



— Omar Khayam 
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Chapter 1 
Introduction 



Since it was first observed, the origin and nature of the intracluster (IC) light, or diffuse 
light, has been a mystery to astronomers. With the advent of more powerful telescopy 
the existence of IC stars became clear (Fig. 1.1) and the question became "Where do 
these stars originate?" These stars could be part of cluster-member galaxies; they could 
have formed in-situ; or they could have been stripped from member galaxies by various 
processes. Dwarf galaxies (DGs) are the most numerous galaxies; they are low mass 
(1O 7 -1O 9 M ) galaxies having an absolute magnitude fainter than Mb ~ —16 to —18 
mag (Grebel 2001; Bothun et al. 1991). Contemporary opinion ascribes the origin of IC 
stars to, for example, the tidal destruction of DGs during the evolution of the cluster. 
Clusters are large groups containing up to thousands of galaxies. DGs can be destroyed 
in various ways during the evolution of the cluster. Understanding the nature of IC 
light is important as, for example, reproduction of the observed proportions (and other 
related parameters) of IC light can be used to refine cosmological models. The goal of 
this research project is to compare the various modes of dwarf galaxy destruction and 
ascertain if dwarf galaxy destruction by tides is sufficient to explain the observed IC 
light. 

As part of a team I have simulated and examined the relationships between six 
possible outcomes of dwarf galaxy destruction during the evolution of isolated clusters. 
Subsequently, with this team, I have analyzed the results to date of an extension of 
the isolated cluster results to a cosmologically significant region of the Universe. The 
analysis, and the simulations I carried out as member of the team in this project use 
known parameters such as the luminosity function (LF), various cluster halo mass 
profiles, and the cosmological parameters to simulate, test, and corroborate the method 
used and the results obtained. 
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Figure 1.1 Trails of gas in a nearby cluster. Trails of gas and stars appear to spring 
from two galaxies (at the left and in the bottom right corner) in the Abell 1367 cluster. 
The trails are probably due to a tidal event 50 million years ago. Image credit: Peppo 
Gavazzi, Universita di Milano Bicocca/MPIA. 
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Typical simulations exploring the evolution of clusters devote much time and pro- 
cessing power to a single simulation. Shifting attention from sheer processing power to 
clever programming, however, can allow for statistical exploration of the relationships, 
in a significant region of the Universe, between different modes of DG destruction in 
clusters. DG destruction has been widely touted as an explanation for the observed IC 
light, which contributes up to 25% to cluster luminosities. The destruction of DGs (by 
tides) would produce IC stars in the IC medium; it is these stars that are the purported 
source of the observed IC light. 

Heuristic modeling of the DG interactions allows prediction of the average DG be- 
havior with an economy of time and processing power. Rather than unnecessarily flex 
processing power "muscles", subgrid physics are implemented in this project and a small 
number of particles (~ 1 000) is used to represent a typical cluster. Results are consis- 
tent with accepted theories of the IC light being largely produced by tidal disruption 
of DGs. 

Alternative theories as to the origin of the IC stars include for example, in situ 
formation, with the notion that the IC stars originate at least partly from DGs that 
have had the outermost material stripped by a process known as 'galaxy harassment' 
(Richstone 1976; Moore 1998), or another phenomenon termed tidal stirring, where 
tidal shocks strip the DGs (Mayer et al. 2001). 

In this research project six possible outcomes are identified from interaction between 
a dwarf galaxy and the rest of the cluster: (1) The galaxy merges with another; (2) The 
galaxy is destroyed by the tidal field of a larger galaxy, but the fragments accrete to that 
larger galaxy; (3) The galaxy is destroyed by tides of a larger galaxy, and the fragments 
are dispersed in the IC medium; (4) The galaxy is destroyed by the tidal field of the 
background halo; (5) The galaxy survives to the present (i.e. it is not destroyed by any 
process); and (6) The galaxy is ejected from the cluster. By analyzing these interactions 
it is shown that the average behavior of DGs is consistent with the destruction of DGs 
as an explanation for the origin of IC stars when one takes into account the various 
results of interaction. That is, IC stars come from a variety of origins, which for the 
most part arise from DG interactions in the cluster. 

The results from this part of the project, obtained from simulations of a single 
idealized cluster are extended in a second, ongoing half of the project that simulates a 
cosmologically significant volume of space, i.e. a volume containing numerous clusters. 

Study of the IC light is timely, as recent discoveries of IC stars in the Virgo cluster 
might result in an opportunity to learn about details of cluster origins (Arnaboldi 2004). 
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Additionally, understanding properties of clusters, such as the luminosity function, 
number density function, and formation epoch is worthwhile as these are important 
tools in testing cosmological models and parameters. 



1.1 Dwarf galaxies 

Galaxies in general can be thought of as the building blocks of the Universe; DGs can 
be thought of as the fundamental units of mass. DGs are the galaxies that have the 
largest dark matter (DM) component (Bothun et al. 1991) and they are the most 
numerous. Nonetheless, most studies conclude the DG contribution to the total mass 
density is smaller than that of giant galaxies, stemming from the untested assumption 
that the mass-to-light ratio (M/L) of galaxies does not increase for smaller and smaller 
galaxies (Bothun et al. 1991). The simulations presented here include only gravity, with 
no other physical process (gas and DM behave the same way under gravity). When 
necessary, the IC light is calculated using a M/L. 

DGs, defined in the previous section as low mass (1O 7 -1O 9 M ) galaxies having an 
absolute magnitude fainter than Mb ~ —16 or —18 mag, have low surface brightness and 
low metallicity (Grebel 2001). In the simulations presented herein, as in the observed 
Universe, DGs are the most numerous galaxies occurring. DGs are thought to comprise 
approximately 80% of the Local Volume galaxy population (D < 10 Mpc), and may 
have a space density ~ 40 times that of bright galaxies in the Universe (Staveley- 
Smith et al. 1992). There exists, however, a deficiency in the number of observed low 
luminosity DGs (discrepancy more than one order of magnitude) as compared with the 
large number of theoretically predicted low mass DM haloes (Trentham & Hodgkin 
2002; Trentham et al. 2006). This is recognized as a problem for cold DM theory; the 
likely solution involves energy feedback from stellar evolution. 

DGs belong typically to one of three main morphological types: (1) dwarf irregulars 
(dlrr), (2) dwarf spheroidals (dSph), and (3) dwarf ellipticals (dE). These types of 
galaxies have similar scaling, but some dwarf irregulars have bluer colors and more 
disturbed morphologies; they are called blue compact dwarfs. There exist other, rare 
galaxies that do not fit into this classification, such as huge low surface brightness giants 
like Malin 1 (Trentham 1998). 

There is evidence that dEs are further subdivided into two classes, compact ones 
and diffuse ones. The only compact dE in the Local Group is M32, which lies very 
close to the Milky Way's giant companion, M31. Compact dE galaxies that are not 
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close companions of any giant galaxy have been found in nearby galaxy clusters. The 
radial surface brightness profiles of these galaxies are reasonably well fitted by a i? 1 / 4 
law, corresponding to their elliptical nature; this is called the de Vaucouleurs law, 
as per Binney & Merrifield (1998), hereafter BM98. By contrast, the radial surface- 
brightness profiles of the other class of dE galaxies, the diffuse ones, are best fitted by 
an exponential e~ R ^ Rd , 

Bothun et al. (1991) mention that there are at least two other profile types. One is 
that some dwarfs have very flat cores and then steep fall-offs; in general these objects 
are at the faint end of the distribution in central surface brightness, and consequently 
they are the most difficult to discover. Another is that some compact dwarfs have 
power-law surface brightness profiles. If sufficiently faint or small, these galaxies are 
difficult to distinguish from star clusters. 

The DGs in the simulations of this project are assumed to be spherical (an idealiza- 
tion). When calculating galaxy-pair interactions, however, a galaxy is herein considered 
a bound virialized system with an internal kinetic energy and a potential energy. These 
energies are included in the total energy of the interacting pair of galaxies (Barai et al. 
2009). 



1.2 Dark matter 

The as-of-yet unobserved dark matter (DM) is a pressureless fluid thought to be respon- 
sible for the growth of cosmological perturbations through gravitational instability; it 
is expected to be more abundant in extensive haloes that stretch to 100 - 200 kpc from 
the center of galaxies (Vergados 2006). 

The single systems with the largest proportion of DM are believed to be faint DGs. 
DGs in general have a correspondingly high ratio of dark-to-luminous mass, with M/L 
as high as that of galaxy groups and poor clusters. DGs are increasingly DM-dominated 
at fainter magnitudes. Among the faintest dwarfs in the Local Group are DGs such as 
Draco and Ursa Minor with M/L »s 100 within their core-fitting radii and total M/L 
that are probably much larger than this (Trentham 1998) 1 . 

It is thought today that of the total matter in the Universe, an overwhelming 95% 
of the mass in galaxies and galaxy clusters is made of an unknown DM component, 

1 The high DM content of these faint DGs suggests a link between their mass functions and the 
power spectrum of primordial fluctuations of small scales (Trentham & Tully 2002). 
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the vast majority of which is non-relativistic (cold). Evidence for the existence of 
DM comes from measurements of the rotation curves of galaxies, measurements of 
gravitational lensing and from the existence of hot gas in clusters (Freese 2006) as 
well as from methods of estimating cluster masses on the basis of the temperature 
of IC gas (Rosati et al. 2002). Evidence also comes from measurements including 
cosmic microwave background temperature fluctuations, redshift surveys, gravitational- 
lensing shear effects, and distribution of absorbing neutral clouds along different lines 
of sight (Pryor 1996). The concordance cosmological model, ACDM (where A is the 
cosmological constant and CDM stands for cold DM), is the simplest model resulting 
in an interpretation consistent with the observed Universe. Figure 1.2 shows the results 
of various simulations of which the ACDM model produces voids and filaments closer 
to observations of the actual Universe. 

The total density parameter is f2 = 0.27. Baryons account for £\ — 0.04. Hence, 
the non-baryonic (exotic) part accounts for 0.23/0.27 = 85% of the mass in the Universe, 
and this component is what is normally referred to as DM. Notice, however that some 
of the baryonic matter is unobserved. This component is called baryonic DM. The 
most likely candidates for dark baryons are known as massive compact halo objects 
(MACHOs), which include faint stars or stellar remnants. The best current limit on 
the MACHO abundance is 20% of the dark halo mass. Interestingly, the appearance of 
baryons varies with scale, occurring mostly in stars below a galaxy mass scale of M* (M* 
is the mass of the Galaxy, the Milky Way) and mostly in hot gas for systems much more 
massive than M* (i.e. galaxy groups and clusters), while for DM there is no difference 
in properties between galaxy, group, and cluster scales (Silk 2007). Further progress on 
DM structure will come as progress is made in understanding when and how galaxies 
formed (Silk 2007). Currently, it is thought that a realistic distribution of DM haloes 
might be produced by models including a process involving feedback augmented by, for 
example, hyper novae (a very energetic supernova thought to result from an extreme 
core-collapse scenario), or a top-heavy initial mass function (the mass dependence of 
the number of stars that form per mass interval per unit volume). Models resulting 
in such realistic distributions can involve for example incorporating star and active 
galactic-nuclei (AGN) feedback; using massive gas outflows can effectively weaken the 
DM gravity in the central cusp of the dense baryonic core, which forms by gas dissipation 
(Silk 2007). 

The presented simulations assume an isotropic distribution of the galaxies inside the 
cluster; this is an idealization, but little is yet known of the internal structure of the 
DM as explained in the foregoing. Furthermore, the positions of the galaxies herein are 
chosen using realistic guidelines consistent with the literature. 
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Figure 1.2 Predictions of various cold DM models. States for three redshifts are 
shown. ACDM is the concordance model; it includes a cosmological constant. The stan- 
dard model (SCDM) does not assume a cosmological constant. The tilted model (tCDM) 
takes into account possible pre-inflationary inhomogeneities. The open model ( OCDM) 
produces results consistent with the observed Universe but does not assume the Universe 
is flat (flatness is more consistent with inflation). Image credit: (GIF) Colberg (1997). 
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1.3 Luminosity functions 



The relative number of galaxies of various Hubble types is usually represented in terms 
of the luminosity function (LF) as shown in Figure 1.3. The galaxy LF is closely 
related to the galaxy mass function, one of the most important parameters in galaxy 
formation models (Trentham & Tully 2002). The LF (loosely based on the Press- 
Schechter formalization for the primordial halo distribution) has for nearly 30 years 
been a standard tool for quantifying the galaxy population (Driver 2004) and is a useful 
description of the galaxy content of any particular environment as it is straightforward 
to measure (Trentham & Tully 2002). Recent studies consistently show that, with some 
exceptions, low surface brightness (LSB) is synonymous with low luminosity (Driver 
2004). Small stellar fraction and very low luminosities make DGs the hardest galaxies 
to detect. It is in fact generally problematic to quantify the low end of the LF as 
illustrated by the lack of points contained under the low luminosity area of the curve 
in Figure 1.4 . 

The LF can be defined as the number density of galaxies per unit luminosity (Tren- 
tham & Tully 2002). A functional form of a general analytic fit to galactic LFs, as 
proposed by Schechter (1976), is: 

$(L) = ($*/L*)(L/L*) a e- L/L *, (1.1) 

or equivalently, in terms of absolute magnitude: 

$(M) = (0.41nl0)$*10 a4(a+1)(M *- M) exp[-10 a4(M *- M) ], (1.2) 

where a and L* (or M*) are free parameters used to obtain best possible fits to available 
data (Carroll & Ostlie 1996). The normalization at the characteristic luminosity L* is 
$*; it is used as such a free parameter in some descriptions of the LF. Subdividing 
the morphological galaxy types finely, shows that different classes of galaxy have very 
different LFs (BM98). The LF consistently provides a good formal fit to the observed 
luminosity distribution. This consistency between the luminosity distribution and the 
LF appears to hold regardless of environment. A value of a = — 1 implies equal numbers 
of galaxies in the magnitude intervals, a more negative (or steep) value implies that 
dwarf systems are more numerous (Driver 2004). 

To explain the shape of the LF in a cold DM Universe, it is necessary to include 
at least a feedback mechanism (beyond the heating resulting from the photo-ionization 
of the pregalactic gas) to flatten the faint end of the LF and to suppress cooling at 
the centers of the massive haloes of groups and clusters. Despite various mechanisms 
having been attempted, however, it is not a simple matter to replicate the shape; for 
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Figure 1.3 The global luminosity function. Shown here is an example LF for a 
nearby volume-limited sample and images of the actual galaxies contributing to the LF. 
Image borrowed from Driver (2004)- 
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Figure 1.4 The lower end of the luminosity function. The yellow box marks the 
boundary of the sample limited at B = 20 mag. This box barely contains volume for very 
low luminosity systems. On the left are seen the distribution of galaxies from the local 
group, Fornax, and the Galaxy's globular clusters. Image borrowed from Driver (2004)- 
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example, to reproduce a sharp cutoff as observed at the bright end of the LF; Benson 
et al. (2003). 

The question concerning the LF at the faint end is an old problem noted by Zwicky 
in a private communication to Parolin et al. (2003): "Where does it [the faint end of 
the LF] end and what kind of objects populate it?" This problem has become easier to 
address with the high resolution imaging and spectra in modern times, and yet "fainter 
than Mr = —16, the LF is only known in a handful of environments" (Trentham & 
Tully 2002). 

Barkhouse et al. (2002), analyze a very large sample of cluster galaxies with the 
same purpose as Parolin et al. (2003) namely to examine the dependency of the faint 
end slope on cluster morphology as well as to examine the gradient, if any, in the 
faint end as one moves towards the outskirts of a cluster. They do not go very deep 
in magnitude (M R < —16), but their results offer positive evidence. There has been 
a suggestion that the faint-end slope of the galaxy LF may be steeper than previous 
studies indicated. This suggestion arises from the existence of very faint DGs not being 
included, such as those having very flat cores and then steep fall-offs (therefore being 
very hard to detect). Including these dwarfs would increase the faint-end slope beyond 
— 1.5. Measurements of tidal radii and carbon star velocity dispersions in Local Group 
dwarf spheroidals indicate that the largest dynamical M/L belong to the least luminous 
galaxies, thus the M/L of LSB DGs could be large. The M/L numbers in the local group 
range from 5-10 for Sculptor to 50-100 for Ursa minor (Bothun et al. 1991). 

There appears to be no clear sign of a lower limit to the LF of irregular galaxies. 
They are, on average, much fainter than the other galaxy types. The numbers of 
irregulars grow as one moves to fainter magnitudes. Their faint nature can be crudely 
described by a Schechter function (parameters Me*(Irr) = —15 + 51ogh and a(Irr) = 
—0.3, where 100 h is the Hubble constant). A fit to the dE galaxies can be given by a 
Schechter function (parameters Me*(dE) = —16 + 51ogh and a(dE) = —1.3). Dwarf 
ellipticals show a LF that is open ended, similarly to the irregular galaxies (BM98). Note 
that in our simulations the Hubble constant is given a value H = 73.2km s^Mpc -1 
consistent with the results of WMAP3. 

Although the field and cluster LFs are broadly consistent (Fig. 1.5), upon comparing 
the respective LFs for the local field of galaxies near the Milky Way and the Virgo 
cluster, it is clear that there is not a universal LF. Rather, each LF depends on the 
environment of the particular sample of galaxies (Fig. 1.6). Similarly observations 
suggest that the galaxy type depends on the environment (Fig. 1.7). Astronomers may 
in the future have the alternative of using information from bivariate distributions as 
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highlighted by Driver (2004), who emphasizes that the global galaxy LF (red line in 
Fig. 1.3) condenses the available information of galaxies (images in the Figure) into three 
crucial numbers: 1) the characteristic luminosity M*; 2) the absolute normalisation $*; 
and 3) the faint-end slope a. Driver (2004) highlights in addition, that although the 
Schechter parameterization is usually a remarkably good fit, it may be that too much 
important information is lost, for instance the sizes and bulge-to-total parameters. 

In their long term research, Parolin et al. (2003) have measured a set of clusters 
in various colors in order to better understand the shape of the LF in clusters and 
to pinpoint the slope of the faint end. Their analysis and data interpretation use 
different algorithms to reduce or eliminate systematic error, robustly corroborating a 
bimodal nature to the cluster LF. They emphasize, however, that in their view, it is the 
estimate of the total luminosity in different optical color and its correlation with the 
X-ray luminosity that will give the opportunity, compared with the X-ray luminosity, 
to better pinpoint some of the characteristics of the cluster and IC medium evolution. 

While there is no universal LF, the luminosity function of cluster galaxies is nearly 
the same from cluster to cluster (Fig. 1.5), generally following the form proposed by 
Schechter (1976). Observing the high luminosity tip of that distribution allows one to 
normalize the overall galaxy LF for the cluster, yielding estimates for both the cluster's 
total optical luminosity and its mass. The number of galaxies in luminosity range dL, 
about L, is proportional to L a e~ L l L * , with a ~ — 1 (Voit 2005). Carlberg et al. (1997) 
showed that the radial mass density p(r) and the radial number density u(r) of galaxies 
are roughly proportional to each other; a finding that is used in this project to assign 
a realistic galaxy distribution for galaxies in the isolated cluster. 

The Schechter LF is used variously in this research project to guide the results. The 
LF is used, for instance, in combination with a Monte-Carlo rejection function to gener- 
ate various initial conditions (i.e. the mass of the cluster halo and the mass of the DGs 
inside the halo). The first part of the simulation uses galaxy particles with distributions 
of masses and luminosities chosen to follow a Schechter luminosity distribution and an 
accepted M/L. The initial position of the galaxies (radial coordinate) is then paired to 
a given galaxy avoiding overlap. Several series include a central giant cD galaxy. In 
the second half of the project presented here, when the simulation has proceeded to an 
advanced redshift, the mass distribution is plotted and the corresponding LF is com- 
pared with observed LFs; this comparison is then used to ensure the obtained results 
are realistic (Barai et al. 2009). 
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Figure 1.5 Various luminosity functions for different environments. Data orig- 
inate from various surveys. Dotted lines are extrapolations. Solid lines are data fits. 
The cluster LFs are normalized to $* = 0.0161. Shown here is a broad consistence 
between the field LF and the cluster LF. Figure borrowed from Driver (2004)- 
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Figure 1.6 Luminosity functions for two samples of galaxies. Image borrowed 
from Carrol & Ostlie (1996). LFs are shown here for samples of galaxies in the vicinity 
of the Milky Way and the Virgo Cluster (top, bottom, respectively) . The total LF in 
either environment is the sum of the individual LFs of each Hubble type. 
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Figure 1.7 Galaxy type frequency in different environments. Fraction of morpho- 
logical galaxy types in different environments (bottom). Local galaxy number density is 
shown as projected surface density (top). This figure was borrowed from Longair (1998). 



1.4 Galaxy clusters and the luminosity function 

Clusters of galaxies are a rich information source about the underlying cosmological 
model; they make possible a number of critical tests. Because structure grows hierar- 
chically, the growth and development of clusters directly traces the process of structure 
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formation in the Universe. Clusters are dynamical systems not yet in equilibrium. Clus- 
ter galaxy velocity dispersion and the virial theorem thus cannot yield an exact mass 
measurement. Instead, other data (i.e. spatial distribution of galaxy velocities) must 
be used to accurately measure the masses of nearby clusters (Voit 2005). 

For distant clusters the data necessary to supplement the mass measurements are 
more difficult to obtain. To fill this gap, simulations of cluster formation can be used for 
calibrating the approximate virial relationship between velocity dispersion and cluster 
mass. Yet another of Zwicky's (1937) important contributions is his idea that cluster 
masses could be measured through gravitational lensing of background galaxies. This 
is one of the primary contemporary methods for measuring cluster mass. The masses 
of clusters now are measured at roughly 10 15 times that of the Sun's. This tremendous 
mass notwithstanding, all of the stars in all of a cluster's galaxies represent only a small 
fraction of a cluster's overall mass. The majority of the mass, in clusters as throughout 
the Universe, is an unobserved component, DM. (It is the DM's gravity that caused 
the baryons to gravitate together into the DM's deep potential wells, there to become 
galaxies.) Clusters contain substantially more mass in the form of hot gas than in 
stars. This hot gas is observable with X-ray and microwave instruments. A cluster's 
total mass turns out to be about seven times the combined baryonic mass in stars and 
hot gas (Voit 2005). 

Models for the X-ray surface brightness of clusters are obtained using the hydro- 
static equilibrium equation (variations include assumptions, i.e. spherical symmetry, 
and isothermality of the gas). Although clusters are dynamical systems that have not 
yet completely equilibrated, they appear to be in approximate hydrostatic equilibrium. 
These models are viable because of the relationship between the gas temperature (ob- 
tained from the X-ray spectrum) and the depth of the cluster's potential well; the deep 
potential wells compress the gas, causing it to emit in the X-ray band (Voit 2005). 

As per Voit (2005), under the assumption of spherical symmetry there is, as a 
function of the radius from the density center, a relation between the gas density p g 
and the temperature T: 



dlnpg rflnT 



d In r d In r 




(1.3) 



T 



The characteristic temperature T^ of a singular isothermal sphere with the same value 
of M(r) jr is given by: 

k B T^(r) = GM(r)^, (1.4) 
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where ke is the Boltzmann constant, fi is the mean molecular weight of the IC gas, 
and m p is the proton mass. The well-known f3 models add the assumption that the gas 
is isothermal (Cavaliere & Fusco-Femiano 1976; Voit 2005). Assuming an isothermal 
velocity distribution and a constant velocity dispersion am, gives: 

dhip g _ pm p d(b _ p dlnp 
dr kT dr dr ' 

where p g is the gas density, and p and T are generalized to the particles responsible for 
the mass distribution. Further, as detailed in Voit (2005): 

(3 = pm v a\ D jYY. (1.6) 

The observed surface brightness profiles of clusters is described quite well by the (3 
models in a given radial range, but give the best fits for rich clusters and a possible trend 
toward lower j3 values in poorer clusters. Most of the observed X-rays must come from 
a relatively small proportion of the IC medium (as the X-ray luminosity integrated over 
radius converges for (3 > 0.5; Voit 2005). The (3 models can underestimate the central 
surface brightness within the core radius r c and overestimate the brightness at s ^> r c 
(Vikhlinin et al. 1999). This is in part because of the assumption of isothermality of 
the IC medium (an idealization), and in part because real cluster potentials differ from 
the King model (Voit 2005). 



There is a generic form for representing density profiles shallower than isothermal at 
small radii and steeper than isothermal at large radii, as corresponding to DM haloes; 
this is: 

p M (r) (xr- p {r + r s ) p - q (1.7) 



where p and q describe inner and outer power-law slopes, correspondingly, and the 
radius r s marks the steepening of the profile (Voit 2005). Some models for the X-ray 
surface brightness include the (3 model, the King model (/3 = 1) (Molinari et al. 1998), 
and the Navarro, Frenk & White model (Navarro et al. 1997), hereafter NFW, as well 
as the Hernquist model (Nipoti et al. 2004). 



The (3 model is given by: 



^gas 



Po 



1 + (r/r c y 



-3/3/2 



where po is the central density. The values of p , r c and (3 are taken from Piffaretti & 
Kaastra (2006a), who give the gas density parameters for 16 nearby clusters. Scaling 
the gas density with the universal ratio of matter (dark + baryonic) to baryons gives the 
halo density: p halo = PDM + Pgas = P ga s^M/^b, where Om and fl^ are the present matter 
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(baryons + dark matter) density parameter and baryon density parameter, respectively. 
The foregoing assumes that the cluster baryon mass fraction follows the cosmic value 
of Q^/Qm', this is expected to be generally true, although precise estimations of cluster 
baryon content have shown deviations from the universal value (White et al. 1993; 
Ettori 2000; Gonzalez et al. 2007). 

The much used NFW model is a special case (£ = 1) of the functional form as in 
Biviano & Girardi (2003): 

P(r)= ( , , , 3 - g - (L9) 

{r/ay (1 + r/a) s 

The NFW model is obtained when the distribution of gas and DM in the background 
halo follow analytical models of the DM density having a functional form as: 

Pdm (0 = Ps 2 , (1.10) 

(r/r s ) (1 + r/r s ) 

where p s is a scaling density, and r s is a scale length (Navarro et al. 1997). The NFW 
profile is often parametrized in terms of a concentration parameter c. The parameters 
p s and r s are then given by 

200c 3 p crit (z) = 25H 2 ( Z y 

Ps " 3[ln(l + c)-c/(l + c)] 7rG[ln(l + c)-c/(l + c)] ' l " ' 

r s = ^, (1.12) 

c 

where p C nt( z ) — 3H 2 (z)/8tiG is the critical density at formation redshift z, and r 2 oo> 
the virial radius, is the radius of a sphere whose mean density is 200p cr ; t (200 times 
the critical density of the Universe at the epoch of formation). After scaling, the halo 
density profile is p halo = p dm ^m/(^m - ^b)- 

The Hernquist model differs only in the exponent of the numerator, as below: 

Pdm ( r ) = ~ 3 • (1.13) 

(r/r s )(l+r/r s ) 



The mass density profile can also be approximated by: 

P(r)= 1A P ° L6 , (1.14) 

\r/r s ) (l + r/r s ) 

with p ~ 5900p c and r s = 370 kpc. While the outer-region slope (a ~ 3) is steeper 
then as given by Moore et al. (1998), both of these profiles agree as to the inner profiles 
(a = 1.4; Lewis et al. 2000). More recent high resolution simulations all show ~ r - L5 
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for the inner profile, turning over to r~ 3 at large radii (Lewis et al. 2000; Moore 1998). 
The above equation 1.14 matches the density profile of NFW at large radii but has a 
steeper central cusp. 

The simulations of this project use the above profiles (/3 and the NFW profiles) to 
model the isolated cluster and subsequently the clusters in the simulated cosmological 
volume. In the first case {j3 model) it is assumed that the DM in the background halo 
follows a density distribution similar to observations of the IC gas. The halo density is 
obtained by scaling the gas density with the universal ratio of matter (dark + baryonic) 
to baryons. The density profile is then integrated to obtain the background halo mass. 
In a second case the halo mass is obtained by integrating over the NFW profile, with 
the assumption that the distribution of gas and DM in the background halo both follow 
the NFW profile. In a future project it might be possible to use a more detailed profile, 
perhaps bimodal, as suggested by the results of several researchers as below. 

Parolin et al. (2003) find consistency between their results and other works in the 
literature confirming the bimodal nature LF of cluster galaxies. Driver et al. (1994), 
studying Abell 963 (z 0.2), in their analysis on a sample that reaches Mr = —16.5, 
similarly find a LF that could be fit by a composition of two Schechter functions, 
or by a Schechter function plus a power law with slope —1 and —1.8 respectively 
(in reasonable agreement with other findings). Supportive results are also found by 
Molinari et al. (1998); and by Barkhouse et al. (2002), who analyze a very large sam- 
ple of cluster galaxies as well as the LF of the first cluster of the sample used by 
Parolin et al. (2003); namely their results corroborate a bimodal LF and a faint end 
slope, a = —1.65 (Parolin et al. 2003). 

It is possible that the LF condenses too much into three simple numbers, thus losing 
valuable information. Driver (2004) stresses that the community should use the legacy 
datasets to supercede the LF with for example two multivariate distributions: (bi- 
variate) the luminosity-surface brightness plane (LSP), and the color-luminosity plane. 
Galaxy bulges and galaxy disks form marginally overlapping but distinct distributions 
in both planes. This indicates two key formation/evolutionary processes: merger and 
accretion. 

The luminosity- surface brightness relation and the luminosity size distribution de- 
note two LSPs. They are each readily transformed into the other, luminosity, size and 
surface brightness being related by: 

A*hlr = M + 2.51og 10 [27rR| LR ] + 36.57, (1.15) 

where fi e is the effective surface brightness, M the absolute magnitude, and Rhlr the 
semimajor axis half-light radius in kpc (Driver 2004). 



Chapter 1. Introduction 



20 



The LSP offer a strong argument for their utility: both planes show that disks 
and bulges form distinct but overlapping distributions thus indicating secular evolution 
of these components (two mechanisms and two timescales). The LF is a useful global 
measurement of cluster properties, but this finding presents a strong motivation to begin 
measuring distinct components independently. The LSP appears to provide a direct 
meeting ground to the numerical simulations. This is a most important result as it is 
the cross-talk between simulations and observations that will ultimately bring about 
the real insights into the processes of galaxy evolution and formation (Driver 2004). 
Driver (2004) further raises three questions as worthy of attention as the community 
engages, as he suggests, in quantifying the evolution of the LSP distributions across 
the entire path length of the Universe: "1) Which wavelength is optimal for structural 
studies of galaxies? 2) How might we push back the boundaries into the dwarf regime? 
and 3) Can we connect structural measurements to the properties of the DM halo?" 

Tests concerning the internal structure of the DM halo are beyond the scope of 
this project, but in other simulations aiming at explaining, for example the origin and 
nature of the cusp problem, the bimodality of cluster haloes, and other features of the 
DM structure, a bivariate distribution such as the LSP might be a useful modeling tool 
in a parallel way as the LF is used herein. 



Chapter 2 

The Isolated Cluster: Method 



With a simple FORTRAN algorithm we simulated (243 simulations) the evolution of an 
isolated cluster of galaxies in order to examine the relationships between the different 
methods of DG destruction. 

Careful selection of initial conditions, parameters, and simplifications was essential 
for the successful implementation of the simulation of the isolated cluster. These sim- 
plifications and assumptions were deemed to be acceptable as they do not significantly 
affect the average outcome of the simulations. 

The assumptions regarding the cluster simulations are as follows: 

•Isolated 

•Completely formed (comprising iV galaxies, represented each by 1 particle, orbiting 
inside a background halo of uncollapsed DM and gas) 

•Experiencing neither mergers with other clusters nor accretion of intracluster matter 

•Comprised of a stationary (i.e. non-responsive to forces) spherically symmetric DM 
halo, the density profile of which does not evolve with time 

•Galaxies will incur one of six possible fates: 

(1) galaxy-galaxy merger 

(2) destruction by larger galaxy tidal field, with fragments accreting to larger galaxy 
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(3) destruction by tides of larger galaxy, with dispersion into the IC medium 

(4) destruction by background halo's tidal field 

(5) galaxy survives to present 

(6) ejection from cluster 

This part of the project simulates isolated, virialized clusters thus the simulations 
are not "cosmological". This not withstanding, the particular choice of cosmological 
model enters the picture twice: in the determination of the radii of galaxies, and in the 
calculation of the elapsed time between the initial and final redshifts of the simulation. 

For the cosmological model the assumptions made are as follows: 

•ACDM model 

•Matter density parameter, Q M = 0.241 

•Baryon density parameter, fib = 0.0416 

•Cosmological constant, fiA = 0.759 

•Hubble constant, H = 73.2 Km s^Mpc" 1 (h = 0.732) 

•Primordial tilt, n s = 0.958, and 

•Cosmic microwave background temperature, T CMB = 2.725K, consistent with the 
results of WMAP3. 

The remainder of this chapter details the author's contribution to the development 
of the simulations of the isolated cluster. The following chapter (Chapter 3) details the 
author's role in the development of a cosmological simulation that extends the results 
of the isolated cluster to a cosmological volume. 
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2.1 Code and tests 

In preparation for the simulations it was necessary to ascertain that the research project 
had not already been carried out elsewhere. I carried out extensive searches in parallel 
with other team members. These searches were carried out using various search en- 
gines such as NASA's Astrophysics Data System (a free tool available on the Internet), 
arXiv.org, and annualreviews.org, among others. Articles deemed to be relevant were 
summarized for subsequent group discussion. This search was continued beyond the 
initial stages of the project. 

While conducting the literature search, the members of the research team examined 
the particle particle (PP) code to be used. The PP code is written in Fortran 77. 
Fortran is used in astrophysics simulations because of its speed with float operations. 
Additionally, legacy code is often written in this language as translating to something 
more modern, such as Fortran 90 would involve much unnecessary work. The standard 
PP code evolves a system of N gravitationally interacting particles using a second-order 
Runge-Kutta algorithm. It was to be used for the simulation of the isolated cluster 
after modifications by the team to include interaction with the background halo, and 
the additions of modules to deal with the subgrid physics. In the resulting, modified, 
algorithm the number of particles iV can vary, as the particles (galaxies) merge, are 
destroyed by tides, or escape the cluster (Barai et al. 2009). 

The system we model in this project is an isolated cluster comprising N galaxies of 
mass rrii, radius Sj, and internal energy U{ (i — 1, . . . , N), orbiting inside a background 
halo of uncollapsed DM and gas. For many simulation runs cluster properties similar to 
well observed clusters (Virgo, Perseus) were used. Galaxies in the code are represented 
using one particle per galaxy. By comparison with contemporary simulations (i.e. 10 6 
particles), this is a very small number. The use of subgrid physics, embodied in various 
subroutines added subsequently, made it possible to use such a relatively small number 
of particles. In turn this permitted the use of a direct method (such as the PP code) 
instead of a more complicated mesh algorithm, for example. 

Before carrying out the actual simulations, it was decided that tests were necessary 
to ensure the simulation would begin on strong footing. These tests were carried out 
for two particles in orbit and for a collapsing uniform sphere. 

Plugging two particles into the PP code yielded a test simulating the two particles 
orbiting each other around a common center of mass. 
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Figure 2.1 Two particle test encounters precession: X-Y plane. Motion of the 
two test particles is shown on the right bottom panel for x-coordinate versus y- coordinate 
positions. Other panels: x-, y-, and z- position versus time. Precession is visible in the 
box showing the trajectory in the X-Y plane. 



The energies and velocities of the particles were as expected, and the particles or- 
bited smoothly around a common center but the orbits precessed when the minimum 
distance of approach between the particles was shorter than a given distance, the soft- 
ening length e (Figs. 2.1 and 2.2). The softening length in this test was necessary to 
prevent an infinite force as these point particles approached each other (in this test 
we had set e = 0.1). This is important in the subsequent calculation of the force each 
particle exerts on other particles during the simulation. 

The results of this test were analyzed jointly. The consensus for this test was that 
the softening length was too large and that this led to particle overlap and the ensuing 
orbital precession visible in the figures. Reducing the softening length (e = 0.0005) also 
corrected a similar problem in the kinetic energy plotted for the particles. 
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Figure 2.2 Motion of the two particles in the X-Z plane. Top left panel: z- 
velocity versus time; top right panel: potential energy versus time; bottom left panel: 
kinetic energy versus time; bottom right panel: total energy versus time. The velocities 
and energies varied in time as expected, with the kinetic energy reaching a peak as the 
potential energy is at a minimum. 



In the second test particles collapsed (in the free-fall time) under gravity, from 
an initial spherical distribution centered at the origin (Figs. 2.3-2.4). Although the 
particles collapsed smoothly, they bounced towards the end of the collapse. At a time 
close to the free-fall time, there was a reversing in the magnitudes of the potential and 
kinetic energies in the related plots (Fig. 2.4). Tuning e eventually produced an ideal 
value of e = 0.001. It was then decided that in the main simulation e would be of order 
of the radius of the smallest galaxies. 
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Figure 2.3 Spherical collapse test; initial positions. This test was done for 503 
particles cubed. The total time is the free-fall time. At t = 0.385 the first particles reach 
the center. At t = 0.3926 the free-fall time is reached. Some particles do not reach the 
center; some particles bounce back. 
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Figure 2.4 Spherical collapse test; energy is conserved in the system. This plot 
shows the evolution of energy of the system of particles under spherical collapse. For 
this test the softening length e = 0.001. The kinetic energy (short-dash line) plus the 
potential energy (long-dash line) add up to zero with energy conservation (middle, solid 
line). 
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2.2 Parameter searches 

Several parameters were required to describe the cluster, these are related to: halo and 
galaxies, softening length, galaxy interactions, and to general simulation components. 
It was decided that cooling flows, accretion to cluster, and galaxy formation would be 
ignored (subsequently a treatment of cluster accretion was implemented). Important 
parameters to find in the literature were: 

• Mass or density cluster profile (from either X-ray luminosity or radial velocities) 

• Number of galaxies to be used 

• Galaxy mass spectrum N(m) (post-search, Schechter LF is used) 

• Relation mass-size for galaxies (observational and theoretical) 

• Initial positions and velocities of galaxies (caveat: massive galaxies near center) 

• Cluster formation epoch 

2.2.1 Density profile 

There was an open question as to which density profile would be used in the simulations, 
but as it was easy to model several by slightly modifying the code, it was decided that 
the /3 and the NFW models for the density profile of the cluster would be used. I had 
also found that the /3 model is widely thought to overestimate the mass, especially at 
large radii. Modeling various density profiles would allow comparison of the various 
results. 



2.2.2 Cluster mass 

The force exerted by the halo is proportional to M(r)/(r 2 + e 2 ) 3 / 2 , where M(r) is the 
mass inside radius s, and e is the softening length of the force. The original plan 
had been to tabulate this quantity as a function of the radius. It became apparent 
that each density profile would need a separate interpolation table for each value of 
the softening length. To reduce the number of interpolation tables it was decided to 
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tabulate the quantity M(r) only, with the division by (r 2 + e 2 ) 3 / 2 being performed 
during the simulation. 

The code to make the tables maketable.f is presented in Appendix B. This code 
was written by Dr. Martel and modified by myself. The code facilitates calculation 
of i.e. the (3, NFW, King, and Hernquist mass profiles. In order to make the tables 
it was necessary to integrate over the function determining the desired profile. In 
maketable.f, Smax represents the distance from the center, out to which it is desired 
to carry out the integration; whereas ds is the incremental step the integration procedure 
will use for the iterations. The number of entries to the table corresponds to npt; this 
in turn gives the level of definition of the produced table. 

In maketable.f, I modified sections of the code that included the algorithm for 
integration. In order to decide on the boundaries (i.e. a distance for Smax, or a limit 
for npt) I conducted several trial runs until the team had decided on the optimal values. 
It was an interesting point in the project to find and define a physical correspondence 
(parameters Smax, ds, npt) between the cosmological object and numerical methods 
for integral calculus. For instance, to correctly produce the interpolation tables, I had 
to define ds = Smax /npt. 

2.2.3 Initial conditions 

It was necessary to choose the initial mass m, radius r, position s, and velocity v of 
each galaxy in order to determine the initial conditions for the simulations. 

The galaxies were assumed to follow a Schechter distribution, such as in equation 1.1. 
From this distribution, a Monte Carlo rejection method was used to select the lumi- 
nosities of the galaxies. The masses were obtained using these luminosities by means 
of an assumed M/L. 

To obtain the slope of the Schechter function of the galaxies in the simulations 
presented herein, it was necessary to plot the relevant luminosities within each run. I 
ran a code that extracted the relevant values and I tabulated the results in a table to 
be presented to the rest of the team. The relevant values were then extracted so as to 
fit the slope of the Schechter function as further detailed in the related article (Barai 
et al. 2009; see Appendix A). 

Once the initial conditions were decided, see Figure 2.5 for an example; and once 
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Figure 2.5 Initial conditions for run A12. Top panel: initial conditions at z — 1. 

The solid circles indicate the virial radii of galaxies. The large circle is the maximum 
distance r = 3Rq = 2.08 Mpc from the cluster center. Lower left panel: same as top 
panel, with symbols rescaled to optical diameter of real galaxies. Bottom right panel- 
enlargement of the central (0.6 x 0.6) Mpc 2 box on lower-left panel. 
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the various modules of the PP code were in order, the last test run became the first 
successful simulation run A. 

The locations of galaxies inside the cluster were selected with the assumption that 
their distribution is statistically isotropic and using observations and conclusions from 
the literature. To prevent overlap, however, and as a certain structure is expected in 
the cluster (i.e. the most massive galaxies probably reside near the cluster center), the 
most massive galaxy was positioned at the center of the cluster. The next seven most 
massive galaxies were placed between radii R and 2R (where R is 3 times the radius 
of the most massive galaxy). The next nineteen most massive galaxies were placed 
between radii 2Rq and 3Ro- The remaining galaxies were placed randomly between 
radii and 3Ro- The radius s of each galaxy is taken to be equal to the virial radius 
r2oo- The initial velocity of each galaxy was obtained by taking the velocity a galaxy 
would have if it were in perfect circular orbit and randomly varying the magnitude 
within a given range. To determine the direction of the velocities of the galaxies, a 
random angle generation technique was followed, similar to that used in generating the 
galaxy positions. 



2.3 Subgrid physics 

The subgrid physics are effected by modules added to the PP code. These modules deal 
with the six possible outcomes of the DG interactions in the simulations here presented. 
These six outcomes are: 

• The galaxy merges with another 

• The galaxy is destroyed by the tidal field of a larger galaxy, but the fragments 
accrete to that larger galaxy 

• The galaxy is destroyed by tides of a larger galaxy, and the fragments are dispersed 
in the IC medium 

• The galaxy is destroyed by the tidal field of the background halo 

• The galaxy survives to the present (i.e. it is not destroyed by any process) 

• The galaxy is ejected from the cluster 

Figure 2.6 shows the geometry of the tides due to other galaxies in the simulation. 



Chapter 2. The Isolated Cluster: Method 



32 



In the code, the subgrid physics subroutines are named: tidesByHalo, encounter, 
tidesByGalaxy, collision, OutsideHalo, and GalxPresent. The module name de- 
scribes the interaction with which it deals. I tested these subroutines, written by Dr. 
Barai, by adding one subroutine at a time, running the code and comparing the ensu- 
ing code outputs. Next the finished module was commented out, the next was added 
and the process was re-iterated. When a given subroutine appeared to be operating 
correctly the next one was added for subsequent testing. The main steps were adding 
mergers, tides, and harassment. 

The amount of data generated in the run-up tests (and hence in the 243 subsequent 
runs) would be prohibitive were it not for the tree-structure directories of UNIX and 
various editing and plotting programs. I carried out tests for the various modules as 
they were written, and organized and managed the storage of this data for plotting and 
comparison. As the various modules progressed I ran the codes to compare the output 
with previous tests. In subsequent runs B I, I varied the active subroutine (with or 
without harassment), cluster properties, and DM profile {(3 or NFW). 



Figure 2.6 Tides caused by other galaxies. Shown is the effect of tides caused by 
a galaxy of mass m 8 on a galaxy of mass rrij and radius Sj. The two largest arrows 
show the gravitational accelerations caused by galaxy i; the two smallest arrows show 
the accelerations caused by galaxy j. 



Once the initial conditions had been chosen, and tests were finished in the previous 
subroutines, the first meaningful simulation run was carried out. Table 2.1 shows the 
results of this series of simulations. 




2.4 First results 
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This series (series A) constituted 16 simulations, using for the background halo a (3- 
profile with (3 = 0.33, a core radius r c = 3 kpc, and a central density p = 8.14 x 10~ 26 gcm 
which is appropriate for a cluster like Virgo (Piffaretti & Kaastra 2006). This series 
does not include galaxy harassment, or a cD galaxy. A total of 10 series were run, 
comprising parameters as summarized in Table 2.1. We analyzed the data over several 
discussions; the resulting explanation is presented in the related article (Barai et al. 
2009; see Appendix A). Some details are excerpted here and summarized in the next 
chapter. 

Table 2.1 shows the run number in column 1, and in column 2 the total mass iVftotai hi 
galaxies at the beginning of the run (in units of 1O 11 M ). Column 3 shows the number 
of galaxies iVtotai- Columns 4 to 8 show the number of galaxies N merge destroyed by 
mergers; the number of galaxies N^ es destroyed by tides caused by a massive galaxy, 
with the fragments dispersed in the IC medium; the number of galaxies N acci destroyed 
by tides caused by a massive galaxy, with the fragments being accreted onto that galaxy; 
the number of galaxies N^f° s destroyed by the tidal field of the background halo; and 
the number of galaxies N e j ect ejected from the cluster, respectively. Column 9 shows 
the fraction by numbers of galaxies / surv that survive to the present. Column 10 lists 
the values of Jm, the mass fraction of galaxies contributing to the IC stars. Column 
11 lists the values of the total luminosity of the cluster coming from IC stars, /ics- 
Column 12 shows the Schechter luminosity function exponent a start . The exponent 
was obtained by performing a numerical fit to the distribution of galaxy masses. The 
masses were determined using a Monte Carlo rejection method therefore the exponent 
can differ slightly from the intended value a = —1.28 (equation 1.1). Averaging over all 
runs yields a start = —1.280 ± 0.020. Column 13 lists a en d, the best fit to the Schechter 
function for the surviving galaxies in the given simulation. 

There is no occurrence of a galaxy being destroyed by tides by the background halo, 
and the number of galaxies ejected is either or 1. There are large variations in the 
other numbers from one run to the next, but some trends are apparent. Typically, 
50% -60% of the galaxies are destroyed. Run A2 is an extreme case, with 78% of the 
galaxies being destroyed. 

The destruction of galaxies by mergers dominated the destruction by tides by more 
than a factor of 2 except for run A7. If the cases of tidal disruption followed by accretion 
are treated as being mergers, then mergers dominated tidal disruption even more. When 
galaxies were destroyed by tides, the dispersion of fragments into the IC medium always 
dominated over the accretion of fragments to the massive galaxy, but the ratio varied 
wildly, from 114 : 7 in run A7 to 66 : 51 in run A8. 
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Column 10 of Table 2.1 lists the values of Jm, the mass fraction of galaxies con- 
tributing to the IC stars, 

/\/f gal -l- M hal ° 

j _ 7 M tides T m tides ^ ]_) 

•A^total — Reject 

where M refers to the mass in galaxies, rather than the number of galaxies. The 
galactic mass contribution to the IC medium consists of galaxies destroyed by the tides 
of another more massive galaxy and by the tides of the background halo (though in this 
series there are no such cases). The galaxy masses are converted to luminosities using a 
mass- dependent M/L from the cold DM cosmological simulations of (Yang et al. 2003). 

The fraction / ICS of the total luminosity of the cluster that comes from IC stars 
(column 11) is calculated as per: 

rgal _|_ rhalo 

f tides ' tides /r> q\ 

/ics — 7 t • ^ > 

Motal — -Reject 



Again, there were large variations; in particular the fraction was very large for 
run A2 and very small for run A5. Averaging over all runs we found 

/ ICS = 0.271 ± 0.095 . (2.3) 



Even though in most cases about half the number of galaxies were destroyed, the 
destroyed galaxies tended to be low-mass galaxies; this explained why / ICS < 1 — / SU rv, 
for all the runs. 

That the galaxies destroyed by mergers and tides or escaping were mostly low-mass 
galaxies leads to a flattening of the galaxy mass distribution function. We computed: 



a end = -1.206 ±0.040. (2.4) 



as the best numerical fit to the Schechter LF. 



Chapter 3 

The Isolated Cluster: Results 



3.1 The series — analysis 



The results of the isolated cluster simulations are being published (Barai et al. 2009, see 
Appendix A). In this article a detailed account is presented of the simulations. Following 
is a description of my contribution to research related to the article, including execution 
of the numerous simulation runs; also following is a summary of the related results of 
these simulations. 

Table 2.2, summarizes the various series of simulations. The first column in Table 
2.2 shows the series label. The remainder of each row details the properties of that 
series. Column 2 gives the number of runs for each series. Column 3 gives the value of 
c^start for each series. Column 4 tells which DM profile was used for each series and which 
cluster that profile resembles. Column 5 gives the value of f3 used, when applicable. 
Column 6 shows the value of the central density. Column 7 gives, for the NFW profile, 
the concentration parameter. Column 8 gives the value of the core radius (the scale 
radius for the NFW profiles). Column 9 tells whether or not each run included a cD 
galaxy. Column 10 shows whether or not each run included harassment. Column 11 
shows, for each run, whether or not there was cluster growth. The series thus represent 
variations of the isolated cluster simulations. I carried out the simulations for series 
A, D, E, F, G, H, I, and J. I also carried out an unpublished (3 profile series with 
«start = 1.34 and two Hernquist profile series with no cD galaxy and a sta rt = 1.28 and 
1.31, respectively. It was jointly decided that the related article would present only the 
NFW and /3 profile related simulations. 



For each simulation run I obtained data that I then tabulated using LaTeX and 
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presented to the rest of the group. Each series thus contains data similar to Table 2.1, 
which presents the results of the first meaningful simulation run. The runs in Table 
2.1 simulate a cluster with no massive cD galaxy at its center and with no galaxy 
harassment. This series uses the /3-profile for the halo, with (3 = 0.33, core radius 
r c = 3kpc, and a central density po — 8-14 x 10 -26 gcm~ 3 . These values are similar to, 
for instance, the Virgo cluster (Piffaretti & Kaastra 2006). 

The data obtained after activating the harassment module are given in series B. 
From the higher tide related destruction rate it was thought that the algorithms used 
credibly modeled galactic harassment while maintaining a Schechter distribution of 
galaxies throughout the simulated evolution of the cluster. 

Variations in the parameters were found to have predictable consequences consistent 
with theory, such as in the increased number of ejections present in series C ensuing 
when the cD galaxy was particularly massive (which in turn increased the radius in 
the criterion for initial positions of the galaxies). Further, using the results of the 
simulations to this point and observations from the literature we calculated the slope 
of the obtained luminosity function in order to average and find the best fits. 

Incorporating a cD galaxy in the simulation (series D) resulted in an increase in 
accretions (i.e. by the cD galaxy) and in a lower fraction of IC stars imparted to the 
IC medium, predictably, as there is less mass available and there are a higher number 
of accretions. 

Parameters such as those of the Perseus cluster (Piffaretti & Kaastra 2006), using 
a (3 profile, resulted in tidal destructions by the background halo (series E); this was 
probably due to the steepness of the cluster density profile in this simulation. This 
conjecture is supported by the specifics of the destructions. A cD galaxy, when present 
(series F), competes for tidal disruptions with the background halo, although mergers 
continue to dominate tides. 

In series G, H, and I simulated a cluster with NFW profile, while varying the slope 
of the LF and the presence of a cD galaxy. In the simulations of NFW clusters, a 
large number of galaxies were destroyed by the tidal field of the cluster halo. The 
number of halo-tide destructions is comparable to, or larger than that of the galaxy- 
tide destructions. This research team theorizes that this is due to the steepness (within 
a given radius) of the NFW halo mass profile, resulting in galaxies near the cluster 
center experiencing a larger tidal field. 

Further in series G, with a lower LF slope and no cD galaxy, it can be seen that the 
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increase in halo-tide destructions leads to mergers exceeding the galaxy-tide destruc- 
tions; combined, however, the tide destructions by galaxy and cluster-halo compare to, 
or exceed mergers. An associated increase in the mass contributed to the IC stars is 
also seen. 

Series H and I (NFW, including harassment, and having a higher slope LF) are in 
keeping with the above observations, although there are more accretions in run I, where 
the cD galaxy is present (similar to series D and F). When the cD galaxy is present 
galaxies approaching the cluster center tend to get destroyed by the tidal field of the 
halo before the cD galaxy can have any effect. This plays a limiting role in the fraction 
of matter imparted to the IC medium. 

In series J the background cluster-halo is allowed to accrete over time, 33% from 
z = 1 to z = 0, as per Wechsler et al. (2002). We find that this leads to a higher 
number of interactions, which in turn decreases the fraction of surviving galaxies. 



3.2 Cluster profile parameter dependence 

With 60 simulations we explored the parameter dependence of the cluster 
(ctstart — —1-36, with harassment, no cD galaxy). We explored the range 0.3 < 
(at fixed r c and po), these being typical values as per the literature. The results 
simulations and the joint analysis are summarized here (for details see Section 
related article of Barai et al. 2009; Appendix A). 

3.2.1 The (3 profile 

An increase in (3 causes a decline of galaxy interactions, hence an increase in the number 
of galaxies surviving to z — 0. With mostly the smaller galaxies merging, the number 
of mergers dominated the galaxy tides, but the merging mass fraction decreased below 
the tidal destruction mass fraction (at (3 < 0.6). As (3 was increased from 0.3 to 0.6, 
however, it appeared that more massive galaxies were merging. At (3 > 0.6 the galaxy 
accretion dominates the mass- fraction. 

The IC stars luminosity comes dominantly from the DGs destroyed by the tides of 
more massive galaxies. Exploring variations of the core radius r c between 10-500 kpc 
(Table 12 in the related article of Barai et al. 2009; see Appendix A) showed that for 
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small r c , mergers outnumber galaxy tides, are comparable to tides at r c rs 350 kpc, and 
that tides dominate at large r c . As r c rises, a decreasing number of galaxies survive up 
to z — 0. This leads to an increase in the mass contributed to the IC medium. When 
r c > 50 kpc, the mass fraction is increasingly dominated by galaxy-tide destruction, 
with the fragments being dispersed into the IC medium. The merged mass fraction, 
however, is significantly smaller than that of tidal destructions as the merging galaxies 
are mostly low mass. As r c increases, galaxy tides increase substantially by number 
and mass, while mergers remain almost constant. 



3.2.2 The Navarro, Frenk, and White profile 

The dependence on the parameters governing the NFW model cluster halo density pro- 
file was explored with 35 simulations (a B tart — —1-31, with harassment, no cD galaxy). 
As per the literature, typical values used for the scale radius are r s = 100-500 kpc and 
for the concentration parameter, c = 3-6. 

We investigated five different values of the scale radius r s within 10 - 300 kpc, keeping 
the concentration fixed at c = 4.5. With increase of r s , a smaller number of galaxies 
survive up to the present, causing greater mass to be transferred to the IC medium. Near 
r s > 100 kpc, the mass fraction is increasingly dominated by halo tidal destructions. 
When r s is small, the number of mergers is higher than that of tide destructions; but at 
r s rs 300 kpc the numbers of mergers, galaxy-tide, and halo-tide destructions become 
comparable. The merged mass fraction, however, is comparable to that of galaxy-tide 
destructions; these are both dominated by halo-tide destruction. 

It seems that as r s increases, lower-mass galaxies are tidally destroyed by other 
galaxies; the galaxy-tide destructions decrease in number but increase in mass. 

Varying the concentration parameter to c = 4 and 6 with a fixed value of the scale 
radius r s = 100 kpc, it was found that the mergers continued to outnumber the tides. 
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3.3 Discussion 



3.3.1 Mergers and tides 



In the (3 model, destruction by mergers dominates destruction by tides; in the NFW 
model they are of comparable importance. 

The relative importance of mergers over tides is arrived at by calculating the average 
(over all runs of each simulation) of the fractions below as per the related article (Barai 
et al. 2009; see Appendix A): 



/■mergers Am er g e + -^accr (I ~\\ 

/destroyed at ' V^'-U 

1 ' destroyed 

/V gal -i- /v hal ° 

rtides tides ' tides fn q\ 

/destroyed — at ' 

i » destroyed 



Where ^destroyed = Emerge + + iV accr + N*f°. 

By contrast, increasing values of the core radius r c of the f3 model or increasing 
the scale radius r s of the NFW model leads to a reduced number of mergers and a 
growing number of tide destructions, for instance, at r c > 400 kpc or r s > 300 kpc tides 
outnumber mergers. 

As per the simulations herein, the destruction of DGs alone cannot explain the 
observed IC light. While most of the galaxies destroyed by tides are dwarfs (i.e. run 
C3), the destruction of a few galaxies of mass M > 1O U M provides more than 60% 
of the IC light. To account for all the IC light, we conclude that some intermediate 
mass or massive galaxies are getting destroyed by the tidal field of the most massive 
galaxy. This is viable as the mass ratios between the most massive galaxy and the 
destroyed high-mass galaxies are of order 3-5. Other possible explanations are: (1) 
clusters contain many times more DGs than predicted by a Schechter distribution, or 
(2) DGs have a much lower M/L than assumed. These are both unlikely. 
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3.3.2 Intracluster stars 

The LF of the observed IC stars in clusters falls well within the predictions of the 
simulations herein. The obtained results indicate that the tidal destruction of galaxies 
(by other galaxies and by the cluster halo) in clusters is sufficient to explain the observed 
fraction of IC light; this is supported by observational studies, which find a significant 
IC light component arising from DGs in clusters. 

Results from the simulations presented indicate that for the f3 and NFW models, 
/ics increases with the mass of the cluster halo, which is consistent with studies in the 
literature (Lin & Mohr 2004; Murante et al. 2004). 



3.4 Conclusion 

Areas to Improve 

Including the dynamic friction caused by the halo would result in galaxies falling into 
the cluster center earlier in the simulation. We believe starting the interactions (and 
destructions) earlier would not have a significant effect on the final results; we believe 
also that the method used is sound, although generation of the initial conditions involves 
some tunable parameters. 

Three possible areas of improvement in methodology are in the treatment of galaxy 
harassment, the simplified approach with which galactic encounters are treated (specif- 
ically the energy dissipation thereby), and our idealization of an isolated cluster in 
equilibrium (although we believe the cluster can be treated as isolated after the epoch 
of its major mergers). 

Results 

The approach used enabled us to perform a large number of simulations and cover a 
large parameter space, while obtaining statistically significant results. In a simulation 
using a large number of particles to represent each galaxy it would be more difficult to 
establish the relative importance of mergers versus tidal disruptions. 



Chapter 3. The Isolated Cluster: Results 



43 



The results and conclusions are as follows. 

• The destruction of DGs by mergers dominates destruction by tides. The simula- 
tions presented in the foregoing show a dependence on the scale/core radius. 

• The destruction of galaxies by the tidal field of other galaxies and by the cluster 
halo is sufficient to account for the observed fraction of IC light in galaxy clusters. 
There is an increase of /ics with the mass of the cluster halo, as also supported 
by the literature. Our estimate of the IC light fraction is fairly robust. 

• In clusters simulating the NFW model, the large number of cluster-halo tidal 
destructions dominates the mass fraction. In discussing this point, we noted a 
possible solution to the cusp crisis of DM halos: the central cuspy region of the 
cluster DM halo could have inelastic encounters with the member galaxies; this 
could have the effect of injecting energy into the halo and erasing the cusp. 

• In the simulations presented above, the presence of a cD galaxy increases oc- 
currences of accretion, decreases cluster-halo tidal destructions, and reduces the 
IC stars luminosity fraction. This is opposite to the trend found in observations 
(where /ics is higher in cD clusters). 

• The destruction of high-mass galaxies (M > 1O 11 M ) is required, as dwarfs alone 
do not contain enough stars to account for the observed IC light, even if they were 
all destroyed. A few high-mass galaxies thus must also be destroyed, although 
the vast majority of galaxies destroyed by tides are DGs. 



Chapter 4 

The Cosmological Volume 



4.1 Method 



The cosmological volume part of the simulation refined the crude assumption of an 
isolated cluster in equilibrium. A standard particle-particle/particle- mesh (P3M) algo- 
rithm (Hockney & Eastwood 1988) was used with comoving coordinates that expand 
along with the Hubble flow and facilitate calculations (see Appendix C). Particles rep- 
resenting DGs were created in high-density regions. These particles were allowed to 
gravitate together and merge to form more massive galaxies that in turn gravitated 
together to form clusters. This required the DG mass to be significantly larger than 
that of the P3M particles. In this way the dynamical range of the simulation allowed us 
to model a cosmological volume of the Universe while resolving DGs. The simulation of 
this second half of the project used the same initial conditions as previous P3M runs. 
While working on the isolated cluster I carried out one such run to generate experi- 
mental initial conditions. As in the first part of the project, one module at a time was 
added to the P3M code in order to conduct tests. 

After a calibrating run the parameters of the simulation runs were modified to meet 
time and computational constraints (Appendix D). It was desired that the galaxies be 
significantly more massive than the DM particles. In the first part of the simulation 
M min = 1 x 10 9 Mq was used (the upper mass range of DGs as defined in section 1.1). 
With this value of M m i n , DGs would be only 7 times more massive than DM particles, 
which was not considered enough. For this second part of the project then, the value 
used was M min = 2 x 10 9 M . This value was 14 times that of the mass of the DM 
particles, which was deemed sufficient. In order to save time, tests were carried out in 
a 40Mpc box, forming 1 or 2 massive clusters with 256 3 particles. 
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4.2 Creating the particles 

In order to create the initial conditions, I set an initial run using a 40 Mpc box and 
256 3 particles with a 512 3 grid. This involved changing settings for volume, number 
of particles, and grid. This is effected by modifying the 'include' files, the zeldo.inc 
and nadia.inc related to the P3M code. The nadia.inc file is changed as per the 
prescription shown in Table 4.1 according to the desired parameters of each run. 

After this I modified the parameters of the run in the par_z . dat file (see Table 4.2). 
This file contains the cosmological parameters governing the simulation as well as the 
run-identifier consisting of three sequential characters that vary according to the given 
run. A working knowledge of UNIX was very helpful to manage and handle large data 
files such as the initial conditions generated. 

To test the algorithm for forming galaxies (and later for the subgrid physics) the 
initial conditions from the 40 Mpc run mentioned above were used, with one time- 
saving modification. The file was modified so as to stop the run when the density 
reached 200p mea n- This was because no galaxy forms between z ~ 25 and z rs 9. The 
density reached 200p mean at z = 8.79. All parameters produced for the galaxy particles 
and the DM particles at this redshift were saved in a file to be used in subsequent tests. 

The specific prescription to create DGs (M = M min ) was discussed a few times: it 
was decided to create galaxies in high-density regions and reduce the mass of neighbor- 
ing DM particles. In this way, the number of particles (DM + DG) would increase, but 
the total mass would remain constant. When a galaxy is created, the code reduces the 
mass of the neighboring DM particles to conserve mass. The velocity of the created 
galaxy is also adjusted in order to conserve momentum. A radius Sj is assigned to each 
galaxy, similarly to what was done in the first part of the project. A check was imple- 
mented to ensure that a galaxy was not created inside the radius of an existing galaxy. 
As we had not kept to a simple spherical background cluster with a particular density 
profile, it was necessary to implement a special way to handle tidal destruction by the 
background matter. More specifics will be available in an upcoming article, concerning 
the cosmological volume part of the simulation (Brito et al. 2009). 

Because in the adopted prescription for the creation of galaxies the DM particles and 
galaxies enter the gravity calculation in exactly the same way, it was decided, in order 
to keep track of the galaxy particles, to expand the arrays handling the DM particles 
and to add the galaxies at the end of these arrays. 
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Table 4.1. Guidelines file nadia. inc 



NP3 


NP = NP3 3 


N 


NC 


CPU on purplchazc 


32 


32768 


64 


23 


less thanl hour 


64 


262144 


128 


47 


a few hours 


128 


2097152 


256 


94 


a few days 


256 


16777216 


512 


189 


a few weeks 


512 


134217728 


1024 


379 


several months 



Table 4.2. Parameter file par z 



value in par_z file 


label in par_z file 


label description (not part of par_z file) 


Abl 


ID code 


run-idcntificr 


0.27 


fio 


fio 


0.04444 


n b 


fib 


0.73 


Ao 


Ao 


71. 


Ho 


Hubble constant 


0.93 


n 


Slope of the primordial power spectrum 


2.725 


T_CMB 


Cosmic microwave background temperature 


24. 


z_init 


initial redshift (z) 


30. 


Box size 


Size of the cosmological volume to be used 


73661519 


Seed 


random seed 


4. 


Cloud Radius 


Used in Particle-Mesh algorithm 1 
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I therefore modified the P3M code as follows: 

1) In order to make the arrays bigger, to accommodate the galaxies, I modified the file 
nadia.inc by adding the line: parameter (ntotmax= 1.2 * np). Note that np is the 
number of DM particles. 

2) I replaced the dimension np of arrays in the code by ntotmax, as necessary. 

3) I then added to the main code the common block: common /galaxies/ ngal. 

4) Next, I searched the code for subroutines that looped over all particles and added 
the common block at the beginning of each subroutine (tantamount to initializing the 
variable). I also replaced np by np+ngal in each such subroutine, for instance, in the 
subroutine pppm. I repeated this process in the main module. One exception to these 
changes is that in the subroutine assmass, the loop was split into two. One loop handles 
the DM and the other handles the galaxies (loops 35 and 135 in the main code). 

5) Additionally, to permit breaking the runs into portions, it was necessary to specify 
in the code the initial number of galaxies (ngal). This number was specified by adding 
to the param.dat file the line: number of galaxies. This line could be read by the 
code after adding, where necessary, the line read(l,*) ngal. 

As a result of these changes a galaxy in the simulations herein can have three possible 
states: passive, active, and inexistent. Galaxies are active when they are created. Then 
one of two things can happen: 

1) An active galaxy can be destroyed by tides: it becomes passive. 

2) Two active galaxies can merge: one becomes inexistent, the other remains active and 
acquires the total mass and momentum of the two galaxies (a tidal destruction followed 
by accretion is treated as a merger). 

In loops like do i = 1, np + ngal or do i = np + 1, np + ngal there were two 
possibilities for keeping track of the galaxies after the interaction: in the P3M code, 
the inexistent galaxies are excluded but the active and passive galaxies are included; 
or, in the subgrid physics, the inexistent and passive galaxies are excluded and only 
the active galaxies are included. In this way the code keeps track of the fragments of 
a galaxy that has been destroyed by tides. At the end of the simulation clusters are 
identified and the passive particles in each cluster are examined to find what proportion 
of IC stars they contain. 
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After further changes to the code I carried out another test run, changing the final 
redshift in the param.dat file. I plotted relevant results of this run including the 
GAA_create file (contains information about the galaxies being created), the fort. 91 
file (contains the total mass at every step), and the GAA_cpu file. Soon, the bulk of the 
galaxy creation module had been tested and the module was completed. 

Next the number of galaxies being produced was adjusted and the modules for 
merger, accretion, and harassment were tuned. 

Galaxy creation issues 

A galaxy is created when the density of DM particles at a given cell exceeds a 
threshold value (p > 200p mean )- The cells in the grid with p > 200p mean all have 
an equal probability of a galaxy being created at that point. If a cell has a density 
significantly larger then 200p mean , it is possible that two galaxies are created there (i.e. 
if after a galaxy has been created there, but the density is still significantly greater than 
200p mcan , a second galaxy can be created there at the next timestep). For conservation 
of mass, once a galaxy has been created at a given point, surrounding DM particles 
(around the center of mass) have their mass lowered by the code. The velocity of 
neighboring DM particles is also adjusted in order to conserve momentum. 

We encountered a problem that the DM particles surrounding cells where galaxies 
were created had their masses reduced into negative numbers. To understand this 
problem, we compared the redshift at which particles started to have negative masses 
with, for example, the time at which mergers occurred, and with the total mass of 
galaxies created by a given redshift. This problem of negative mass particles persisted 
as successive modules were activated. Various alternatives were explored in successive 
runs, for example, varying bias. The quantity bias was introduced in the P3M code 
to control the probability that a galaxy particle is created at a given cell. Peak cells 
are ascribed a certain probability and a galaxy is created, or not, depending on that 
probability. The quantity bias has a factor adjustable according to time elapsed in 
the simulation. Initially bias was set at 6.0 x 10 3 , according to the elapsed time in 
computational units (1 time unit = 7.7646 x 10 4 Myr). This value of bias makes the 
cell probability ~ 0.5 at z 9. 

The above problem was explained by the fact that galaxies were being created 
repeatedly at, or near the same points. When this occurred, the neighboring particles 
had their mass decreased successively, causing them to have negative masses at times. 
On subsequent test runs, in order to determine how to constrain the value of bias, 
the total mass in the galaxies created was plotted versus the redshift from the file 
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GAA_gal _Num_Mas s . 

The problem was finally corrected by changing the probability weighting. In the 
code this was done by replacing the line weight (nbnb) = exp(-0.5*dr2/w2) by the 
line weight (nbnb) = pmass (m) *exp(-0 . 5*dr2/w2) . 

During the efforts to solve the negative mass problem, and in order to select the 
value to be used for f2 ga i, we experimented with various bias factors by running several 
test runs up to z = 6. I plotted the total mass in galaxies versus redshift in order to 
determine the appropriate bias for creating a given number of galaxies. To obtain these 
plots, I used an IDL code that plotted the desired quantities from the GAA_gal JlumJlass 
file related to the P3M code. 

For these various runs we compared the amount of mass being created in the galaxies 
by plotting versus redshift for different values of bias. A run was executed, commenting 
the bias factor out, to determine the maximum mass that could be produced (Fig. 4.1). 
The value of bias was then adjusted such that the mass ending up in galaxies would 
total f2 ga i ~ 0.03. This test was repeated with this value of bias. The target was 
to have about 10% of the final mass in galaxies, and the rest in P3M particles. In 
order to monitor the test runs it was frequently necessary to make plots to compare 
characteristics of the tests and the simulations; for example Figure 4.1, with a bias value 
of 0, served to gauge if enough matter was being created, whereas Figure 4.2 served to 
monitor the distribution of the matter in the halo. 



4.3 Test runs for subgrids 

Testing the subgrid physics in the cosmological volume part of the simulation followed 
a course analogous to that in the isolated cluster. The subroutines were added one 
at the time and tested until it was deemed each was working as desired. We resolved 
several problems in the course of testing. 

• The subroutine accrete checks, for each galaxy, if there are DM particles inside 
the radius of the galaxy. Particles within that radius are accreted. The mass, 
center-of-mass position, and center-of-mass velocity are conserved. The old DM 
particle still exists, but its mass set to zero. When the module accrete was 
incorporated into the code, I carried out runs and made plots for comparison 
with previous tests. 
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Figure 4.1 Total mass at a given redshift. This represents part of the experimenta- 
tion carried out to determine the desired value of bias. In this figure bias = 0. 
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Figure 4.2 Halo mass versus geometrical factor. This plot was generated to de- 
termine if the geometrical factor used could reproduce a reasonable halo mass. The 
linking-length (0.171 in this plot), is important when determining the particles that 
belong to the halo. 
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The test run including the accrete subroutine was excessively slow, i.e. it took 
327 steps to reach z = 6 (29 hours) compared with previous runs without the 
accrete subroutine (92 steps to reach z — 6, duration < 18 hours). This problem 
was solved by changing the timestep subroutine to exclude from the calculation 
of the timestep particles whose mass had been set to 0. This increased the length 
of the timesteps and reduced the run duration. 

• Analogous problems were encountered in the subroutine creategal: it allowed 
galaxies that had been destroyed to continue having interactions with other galax- 
ies (their mass had been diminished to zero but the particles had been kept to 
keep track of the IC medium). This problem was fixed by excluding the massless 
particles from those loops after the merger had been processed. 

• The same problem as above was encountered in the tidal destruction subroutine. 

• There appeared to be runaway increase in sizes of galaxies. At times galaxies 
attained almost l/10th of the total mass inside the box. Subsequent multiple 
harassments at the same timestep caused the galaxy to experience an uncontrolled 
increase in size, beyond the size of the box. A size cutoff was stipulated for the 
galaxies in the simulation by incorporating a cutoff size for the galaxies in the 
parameter sizMAXcgs initialized in the subroutine ComovToPhysical. 

As a side effect of excluding the massless particles, there were fewer interactions 
in the test runs. In contrast to the isolated cluster, where galaxies grew by mergers 
only, in this, the cosmological volume part of the simulation, galaxies could grow by 
accretion. This resulted in a lower number of mergers. A test run was conducted after 
the exclusion of the massless particles. The test run created 9682 galaxies with a total 
mass of 0.3079 (including those tidally disrupted and those going to IC light). 

A histogram of the mass distribution of galaxies in the run including accretion 
showed that the total mass in the galaxies followed the prescription used herein but 
that the mass distribution did not follow the Schechter distribution (therefore was not 
correct according to observations). Indeed, in the obtained distribution there were very 
few DGs. This problem was addressed by fixing a threshold lower than 200p mca n in the 
galaxy formation criterion. 

The test runs had begun, initially, in a 40Mpc box (M box m 2.4 x 10 15 M ). While 
these runs yielded a reasonable cluster mass, it was necessary to redo several test runs, 
because of the change in galaxy formation criterion to a threshold lower than 200p mean 
(as a fix for one of the bugs). It was also considered that in view of time pressures 
future test runs would be effected in a 20Mpc box with 128 3 particles. 
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The new test runs thus used a 20Mpc box (total mass Mb ox = 2.95 x 1O 14 M ). A 
problem ensued: the simulation produced only a very small mass of IC stars. The vast 
majority of the tidal disruptions led to accretions, with little material going to the IC 
medium. It was thought that, as the IC light fraction increased with increasing cluster 
halo mass, this problem would disappear when simulating a much bigger box. Another 
possible explanation was that, with the current small volume used in the simulation, the 
galaxies never encountered the gravitational potential of any massive cluster (since the 
clusters forming in a 20 Mpc box are much less massive) and their resulting velocities 
were thus smaller. 

It was also thought that some galaxies had simply been created very close to each 
other, their gravitational potential energy thus being higher than their kinetic energy, 
leading to quick accretions. It was considered that if the merger condition were loosened 
by checking for mergers whenever = Si + Sj, some of the accretions might be counted 
as mergers, thus solving the problem. 

A test run was effected at this point, with the parameters modified as detailed above 
and including the subgrid physics. The time for this full test run was ~ two weeks. 

As the tests progressed the mass distribution plots revealed that the simulation was 
beginning to produce realistic results. Figure 4.3 was used to monitor the shape of the 
corresponding luminosity function, and Figure 4.4 was used as another way to monitor 
the distribution of matter and voids. 



4.4 Conclusions 



The cosmological part of the simulation presented herein is at a very exciting point at 
present. The modules for subgrid physics are working in a realistic way. The latest 
test run produces a Schechter-like distribution (Fig. 4.3). This test run includes the 
galaxy creation module as well as the mergers module. A slice of the distribution of 
particles in the 20 Mpc volume is shown in Figure 4.4. This shows a slice of thickness 
corresponding to 4 Mpc of a test run in a 20 Mpc box. There is an ongoing problem: 
the remaining module, for tidal disruption by the background halo, does not work as 
desired in the test runs. As mentioned in the foregoing, there is currently not enough 
matter being imparted to the IC medium. 



Chapter 4. The Cosmological Volume 



54 



1 



I 2 
o 



- 



-» ' 1 1 J 1 ' • T 1 1 1 1 i J T 1 1 1 J ' r * r 



-2 




_i i i i i ' ■ 



Log (Man / 10" 



2 



3 




Log (Luminosity / 10" L*) 

Figure 4.3 Schechter function reproduced. Top panel: mass distribution of the 
galaxies atz = 0. Bottom panel: luminosity distribution of the galaxies at z = 0; (dotted 
curve: Schechter luminosity function with a s tart = —1-28; dashed curve: Schechter 
luminosity function with a stait = —1.10). 
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Figure 4.4 Cosmological volume slice. Final state at z = of the simulation that 
produced the luminosity function shown in Figure 4-3- This run used a 20 Mpc box; the 
plot shows the location of the P3M particles in a slice of thickness 4 Mpc. 
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The reason for this could be the small size of the test box used and the correspondingly 
small amount of matter in this box. It is thought that there is a number of accretions 
that would be counted as mergers, were the galaxies moving in a more massive system 
(the merger criterion was satisfied too easily in the reduced volume). Were the galaxies 
experiencing interactions, i.e. encounters, or background halo tides, they would have 
larger kinetic energy (hence larger velocities). Tests while correcting this problem are 
expected to be time-consuming beyond the time allotted for a master's degree research 
project, hence I obtained permission to present the results to-date. Further progress will 
be presented in an article (Brito et al. 2009). The Schechter distribution reproduced 
from tests until now, combined with the success of the merging module, and previous 
encouraging results from the isolated cluster (agreement with accepted theory) make us 
optimistic that this problem will shortly be resolved. We will then be able to launch the 
main 80Mpc simulation, expected to take ~1 month to complete. Ideally, the analysis 
of that simulation will show that the results from the isolated cluster carry over to an 
actual cosmological volume where non-spherical clusters have evolved under gravity. 

The next step is to 'turn on the tides' (include the module that keeps track of tidal 
forces by the halo) and carry out test runs that will tell us if we are ready to expand 
the test volume to the full 80 Mpc cosmological volume. These test runs were originally 
to be carried out in a 20 Mpc volume, but will instead be carried out in a volume large 
enough that the total mass in the box will be at least comparable to the mass of a 
cluster. In order for the numbers obtained to be significant statistically, the final run 
must be carried out in a volume large enough to form several clusters (~100Mpc). 

The results for the isolated cluster simulation (part 1 of this project) were in agree- 
ment with the fraction of IC stars currently thought to be due to tidal disruption. 
If this cosmological volume simulation (part 2 of this project) continues to reproduce 
numbers comparable to the accepted theory, and if the distributions of M/L continue to 
reproduce a Schechter function it will be strong support for the validity and significance 
of the method used. 
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Article 1: The Fate of Dwarf 

Galaxies in Clusters and 

the Origin of Intracluster Stars 

This article was accepted in March, 2009, for publication in the Journal of Astrophysics 
and Astronomy, Volume 30, Number 1, pp. 1-36. It is also available at arXiv.org: 
arXiv:0707.1533 



The text of the article was written by Dr. Martel, while the subroutines added to 
the code that modeled the simulations were written by Dr. Barai. My contribution is 
as detailed in Chapter 2 of this thesis. 



Appendix B 



Interpolation tables 



This appendix constitutes the latest version of the procedure that calculates the in- 
terpolation tables from which the mass of the cluster halo is obtained for use in the 
simulations. This is a routine written in FORTRAN; the asterisk at the beginning of 
lines denote comments. 



* This program builds interpolation tables for the density 

* profile of haloes. 
* 

* To calculate values for a different density profile simply 

* 'comment out' the current profile and remove the comment marker 

* for the desired profile. The alternative profiles are located towards the end 

* of the user-defined function at the end of this procedure (maketable) . 
* 

* This procedure comprises three parts. The first part creates the table 

* and calculates table entries using an integral calculated in the second part . 

* The second part uses a user-defined function created in the third part. 
* 

program maketable 

* 

* Code written by Hugo Mart el, modified by William Brito, 2007. 
* 

implicit double precision (a-h,o-z) 
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external g 

parameter (pi=3 . 1415926535) 

npt=10000 
smax=20 . 

ds=smax/npt 

* Open output file 

open(unit=l ,f ile='prof ilel . out ' , status=' unknown' ) 
1001 f ormat(5x, lpel3.6) 
a=0 

write (1, 1001) smax 
do i=0,npt 

s=f loat (i) *ds 

call integ(a, s ,g, area) 

f =4. *pi*area 

write(l, 1001) f 

enddo 

close (unit=l) 

stop 

end 

*========================================================= 

subroutine integ(a,b,f ,area) 

* 

* This subroutine computes integrals using Simpson's rule. 

: \ i : ' double precision (a-h,o-z) 

data tol /l.e-05/ 

if(a.eq.b) then 
area=0 . 
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return 

endif 

areaold=0 . 
sumend=F (a) +F (b) 
sumeven=0 . 
sumodd=0 . 

do n=l,25 

i=2.**n 

h=(b-a)/dfloat(i) 
sumodd=0 . 
do j=l,i-l,2 
c=a+j*h 

sumodd=sumodd+F ( c ) 

enddo 

area= (sumend+4 . *sumodd+2 . *sumeven) *h/3 . 

if (dabs ( (area-areaold) /area) . It . tol) return 

sumeven=sumeven+sumodd 

areaold=area 

enddo 

write (6, 1000) 
1000 format(/5x, 'Error, no convergence.') 
stop 

end 



f unc : i or: G(x) 

* 

* This function calculates the integrand, and passes it to 

* the procedure integ to integrate. 
* 

* the integrand consists of the product of a density function (which is 

* a function of x) and x squared, so as to integrate over a volume 

* and obtain the mass . 
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* G(x)=f (x)*x*x 
* 

* 

implicit double precision (a-z) 

if(x.eq.O.) then 
g=0 

else 

* NFW 

* To calculate values for a different density profile simply 

* 'comment out' the current profile and remove 

* the comment marker for the desired profile. 
* 

* g=x*x/(x*(l .+x)**2) 

* Hernquist 
* 

* g=x*x/(x*(l .+x)**3) 
* 

* xi model (using the best fit values xi=1.4 and a=.33*r200 

* from Biviano and Gerardias. There is another value of xi also 

* given there, check it out now!) 

* Note xi=l corresponds to NFW profile. 
* 

* g= ( (r/a) **-l . 4) * ( (1+r/a) ** ( 1 . 4-3) ) 
* 

* King profile (Beta=l) 
* 

g=x*x/(l .+x**2)**(3/2) 

* 

* The beta model (in Biviano and Gerardi this is credited to 

* Cavaliere and Fusco-Femiano (1978) but 

* it is worthwhile to check the original reference!) using 

* a best-fit value of .8 from Biviano and Gerardi as well as several of 

* Piffaretti and Kaastra's best fit values for clusters listed on their Table 2. 
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endif 

return 

end 
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Coordinate conversion 



The particle-particle/particle- mesh (P3M) code used in this project uses comoving co- 
ordinates, not physical coordinates; the computational volume therefore expands with 
time. Although this means that the simulation box is expanding, the galaxies them- 
selves are not expanding. The box expands along with the Universe, so that the mass 
inside the box remains constant. The box gets larger and larger relative to the galaxies. 
If coordinates are used in which the box had a fixed size of 1, then in those coordinates 
the galaxies would shrink with time. In physical units the box size and the initial size 
of galaxies scale like 1/(1 + z), so in computational units, they are both constant. The 
box expands with time, but the galaxies do not. In computational units, the box is 
fixed and the galaxies shrink. 

The P3M code used utilizes a special form of comoving coordinates, denoted by a 
tilde, called supercomoving variables (Martel & Shapiro 1998). To convert the peculiar 
velocities to actual velocities it is necessary to convert the coordinates before calculating 
quantities such as the kinetic energy. 

In supercomoving variables there is a precise normalization for the scale factor. The 
normalization depends on the particular cosmological model. For the models considered 
here (where f2 is the density parameter, a is the scale factor, Ao is the cosmological 
constant, and the subscript '0' denotes the values at the present time) the solution of 
the Friedmann equation and the present value of the scale factor are, as per Martel and 
Matzner (1999): 

1. Einstein de Sitter model: (Q Q = 1, A = 0): 
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a — i , a = 1. 



(C.l) 



2. Open models (Q < 1, A = 0): 



t 2 - 1 



a 



l-^o 



(C.2) 



3. Flat models with non-zero cosmological constant (fi + A = 1) : 



t = 



1 r a 



dy 



2Jl y 3 /2(l + y 3) 1/2 



, a 



Q . 



1/3 



(C.3) 



The solutions for a(t) do not depend explicitly upon the cosmological parameters 
(these are absorbed in the definition of ao). Therefore, for all the models included in 
the database (the database as in Martel & Matzner 1999), there are only three different 
solutions of the Friedmann equation. 

The positions s and velocities u in physical units are given by: 



r = 



1 + z' 



(C.4) 



u 



J box 



H{z)~r nl /2 H (l + z)v 



l + z 



+ 



(C.5) 



where H(z) is the Hubble constant as a function of redshift (H denoting its value at 
the present epoch). The size of the box being considered is L box and v is the peculiar 
velocity. 



Appendix D 

Box Volume and Particle Size 



In the run used as guide to determine the timeframe for the project, the number of 
particles used in the cosmological volume (120 Mpc box) was 128 3 . The mass per particle 
this implied 1 was ~ 3 x 10 10 M — approximately 30 times the minimum galaxy mass 
used in the isolated cluster simulations (M min = 1 x 10 9 M ). 

To avoid carrying out a simulation that would last several months we could not 
simply increase the number of particles arbitrarily. At that time a final run had been 
envisioned using 512 3 particles. This would reduce the mass to 4.86 x 10 8 M (smaller 
than M min by a factor of 2). 

Using an 80 Mpc box the particle mass would be further reduced to 1.44 x 1O 8 M 
(about 1/7 Mm in), which would allow increasing M min to 2 x 10 9 M (a factor of 2). 
This M m j n 14 times larger than the particle mass used in the isolated cluster simulations 
was deemed to be a good choice. The estimated time for such a run to complete was 
~1 month. With these modifications it was decided to 'tune' the code running tests 
with a box size of 40 Mpc; the 1 or 2 massive clusters that would be produced were 
deemed adequate. In tests, 256 3 particles would be used to reduce running time. For 
this volume, the mass scaling yielded 256 3 particles. These tests would require only a 
few days and would be carried out as the code advanced. For example, a test run was to 

1 In computational units the mass of the box is 1, the dimension of the box is 1 x 1 x 1, and thus 
the mean density is 1. Since the mass is 1, for NP equal-mass particles, the mass per particle is 1/NP 
(as in the zeldo . inc code). In physical units the box is ibox X -^box x -^box- The volume is thus L^ ox . 
The mass inside the box is the volume times the mean density inside the box. (The box represents a 
cosmological volume so its mean density is equal to the mean density of the Universe yO me an-) Thus: 

Mbox — Pme&n 

x -^box = x Pcrit x L^ ox = 3 x H§ x fi x L^ ox x l/87rG. Note that as usual G is the 
gravitational constant, Hq the Hubble constant, and f2o the density parameter. Dividing the mass of 
the box by 128 3 gives the mass per particle. 
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be carried out once the addition was effected of the code that accounts for tidal force by 
the background matter (the cluster was no longer restricted to being spherical). These 
tests would serve, for example, to ensure that there were enough galaxies being formed 
in the simulation and to ascertain that the luminosities of galaxies created, actually 
followed a Schechter distribution (Fig. 4.3). 
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